################################################################################
###############08.Table_2.R#####################################################
################################################################################

pols <- readRDS('pols_bjps')

#####Fixed Effect Model#####

library(plm)
library(pcse)
library(lmtest)
library(dplyr)

#####What is a realistic effect?#####
countries <- c('Austria', 'Belgium', 'Czech Republic', 'Denmark', 'Estonia',
               'Finland', 'France', 'Germany', 'Great Britain', 'Greece', 
               'Hungary', 'Ireland', 'Italy', 'Netherlands',
               'Poland', 'Portugal',  'Slovakia', 'Slovenia', 'Spain',
               'Sweden')


pols_realist <- pols %>%
  dplyr::select(country.name, gini_disp, scale_ggini_inc1) %>%
  filter(country.name %in% countries)%>%
  na.omit %>%
  group_by(country.name) %>%
  summarize(maxi = max(gini_disp, na.rm = T),
            mini = min(gini_disp, na.rm = T)) %>%
  mutate(ran = maxi - mini)

summary(pols_realist$ran) # Average range of increase of the Gini coefficient of 1.595

#####Make table of predicted diffrences#####

pols$gini_real <- pols$gini_disp/1.595

library(plm)

real <- plm(scale_econsdw ~ gini_real + 
              scale_ggini_inc1 + 
              scale_socioeconomicsal +
              gini_real : scale_ggini_inc1 + 
              gini_real:scale_socioeconomicsal +
              scale_ggini_inc1: scale_socioeconomicsal +
              scale_ggini_inc1:scale_socioeconomicsal : gini_real +
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")


library(MASS)
library(gridExtra)

#####Data for Predictive Effects Table#####

#simulate 1000 draws from a multivariate random normal distribution based on the model#
test <- as.data.frame(mvrnorm(1000, mu =coef(real), Sigma =vcovHC(real)))

#vector of sorting from observed min to observed max#
nums <- seq.int(from = -1, to = 1, by = 1)

#create data frame for storing output#
marg <- data.frame(x = nums)

#Salience of economic issues at -1
set.seed(1004)
#Simulate marginal effect based on draws from the model#
for(i in 1:3){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*-1 + test[j,10]*marg[i,1]*-1
  }
}

#Vectors to be filled#
sort <- NA
pe <- NA
lb <- NA
ub <- NA


#Make CI based on empirical quantiles

for(i in 1:3){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

#Combine into a singe dataset

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

print(marg.plot.dat)



#Salience of economic issues set at 0 #

for(i in 1:3){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*0 + test[j,10]*marg[i,1]*0
  }
}

sort <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:3){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

print(marg.plot.dat)


#Salience at 1#

for(i in 1:3){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*1 + test[j,10]*marg[i,1]*1
  }
}

sort <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:3){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

print(marg.plot.dat)

#####These values are then combined in to Table #####